Hysteresis multicycles in nanomagnet arrays 
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We predict two new physical effects in arrays of single-domain nanomagnets by performing sim- 
ulations using a realistic model Hamiltonian and physical parameters. First, we find hysteretic 
multicycles for such nanomagnets. The simulation uses continuous spin dynamics through the 
Landau-Lifshitz-Gilbert (LLG) equation. In some regions of parameter space, the probability of 
finding a multicycle is as high as ~ 0.6. We find that systems with larger and more anisotropic 
nanomagnets tend to display more multicycles. Our results also demonstrate the importance of 
disorder and frustration for multicycle behavior. Second, we show that there is a fundamental dif- 
ference between the more realistic vector LLG equation and scalar models of hysteresis, such as 
Ising models. In the latter case spin and external field inversion symmetry is obeyed, but in the 
former it is destroyed by the dynamics, with important experimental implications. 

PACS numbers: 61.46.-l-w, 75.60.Ej, 85.70.-w 



I. INTRODUCTION 



Hysteresis in magnets [1|, |3 is a paradigm for all his- 
tory dependent behavior in nature 0, 0, [J- In addition, 
hysteresis is the cornerstone of the magnetic data stor- 
age industry, and of great technological importance 0. 
Understanding the full possibilities of magnetic hystere- 
sis is thus important both for fundamental insights and 
practical implications. 

Recently, we have shown that spin glasses can ex- 
hibit stable "multicycle" hysteresis loops, in that when 
the external magnetic field is cycled adiabatically, the 
magnetization returns to itself after m > 1 cycles of the 
magnetic field. This behavior should be experimentally 
observable in spin glass nanoparticles at low temperature. 
Thus simple one-cycle hysteresis loops, although ubiqui- 
tous and generally assumed to be universal, are only part 
of a much richer phenomenon. 

In our previous work l'^, we used the standard 
Edwards- Anderson spin glass Hamiltonian 's\ , with Ising 
spins and nearest neighbor interactions in three dimen- 
sions. Zero temperature dynamics or Monte Carlo dy- 
namics at low temperature were used. Starting from 
saturation, the magnetic field was lowered adiabatically, 
and then cycled repeatedly over a suitably chosen range. 
When the system reached steady state, the order of the 
multicycle (as defined in the previous paragraph) was 
measured. Whether a multicycle is present, and if so, its 
order, depended on the realization of randomness, vary- 
ing from one system to another. 

There were two weaknesses in the previous work. First, 
the systems considered were susceptible to thermal noise 
meaning that very low temperatures were needed to pre- 
vent fluctuations from destroying the periodicity. Sec- 
ond, the model considered was an Ising model. Although 
this has well defined equilibrium statistics, it is not clear 
that it adequately represents the dynamics. 

In this paper, we overcome the problems of the pre- 
vious paragraph by examining a different system as a 
candidate for multicycle hysteresis behavior: an array of 



magnetic nanoparticles. With current technology, it is 
standard to fabricate such magnetic arrays according to 
specification, and a wide variety has been studied exper- 
imentally and these have important applications for 
magnetic storage technology, so called "patterned me- 
dia" Foi' our purposes, these systems offer many 
experimental advantages over the spin glass nanoparti- 
cles that we considered earlier. One can build "designer 
arrays" with optimized parameters that maximize multi- 
cycles and select them in a predictable manner. Spatially 
resolved measurements are possible, unlike for the spin 
glass case, making experimental observations of multicy- 
cles more practical. The relatively large size and shape 
of these also greatly reduces the effects of thermal noise, 
(which is one important reason that they are useful in 
disk drive recording) . Finally, the ability to address small 
regions of the array instead of only applying a uniform 
external magnetic field, and the sensitivity of multicycle 
phenomena to relatively small changes in system param- 
eters, may open up the possibility of using these systems 
for computation. 

Through numerical simulations, we demonstrate that 
multicycles can be seen in an array of pillars made of 
ferromagnetic material, coupled to each other through 
dipolar forces, arranged in square and triangular lattices. 
The external magnetic field is applied perpendicular to 
the lattice, i.e. parallel to the axes of the pillars; to be 
specific, we will refer to this as the vertical or z-direction. 
The existence of multicycles is robust, persisting over a 
large range of system parameters. 

We also examine the importance of frustration and dis- 
order in achieving multicycles, a question that has been 
of great interest in spin glass research Q. We find that 
when a square lattice is used instead of a triangular one, 
multicycles are not as likely to be seen. As discussed later 
in this paper, the magnetization in the pillars prefers to 
be approximately vertical, so that the dipolar coupling 
between the pillars is antiferromagnetic. For a square 
lattice, the dipolar forces between nearest neighbors are 
not frustrated, unlike the case for a triangular lattice. 
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(There is still some frustration because of further neigh- 
bor interactions.) This result emphasizes the importance 
of frustration. 

For the spin glass system considered earlier, disor- 
der was explicitly present through the random bond 
strengths. For the nanoparticle array, although the bond 
strengths are not random (unless the spacings between 
the pillars are varied), there is crystal anisotropy, aris- 
ing from the fact that the magnetization prefers to align 
itself in a specific direction relative to the crystal axes. 
Because of the way in which the pillars are grown, the ori- 
entation of the crystal axes is different in each pillar, and 
random. Even if the crystal anisotropy energy is small 
compared to the dipolar (and other) energies, we find 
that it is sufficient to cause multicycles in an otherwise 
regular triangular lattice. The order of the multicycle for 
a specific sample depends on the orientation of the crys- 
talline axes in its pillars. Although in this paper we only 
consider the case of randomly oriented crystal axes, if the 
pillars could be grown with the orientations specified, it 
would be possible to make arrays whose hysteresis loop 
are multicycles of desired order. 

Regardless of the source, it is desirable to have some 
inhomogeneity in the model to see multicycles. Other- 
wise, when the external magnetic field is reduced from 
saturation, the magnetic moments would be stuck in un- 
stable equilibrium; in practice, thermal noise would push 
the system out of this unstable point. Although this 
effect would not be so extreme subsequently when the 
external field is cycled over a range that does not reach 
saturation, it is clear that to reduce the sensitivity to 
thermal noise, and thereby increase definite multicycle 
periodicity, explicit quenched disorder is desirable. 

The model used for the dynamics in our numerical sim- 
ulations is discussed in detail in the next section; the 
magnetic moment of each pillar is treated as a single 
'Heisenberg spin', i.e. with its orientation as a continu- 
ous variable, with continuous time dynamics. This is in 
contrast to the earlier spin glass work, with Ising spins 
and discrete (event driven) dynamics. Even though the 
dynamics are continuous, as discussed in the next sec- 
tion there is a shape anisotropy energy for tall pillars 
that causes the magnetization to be nearly vertical, and 
jump from up to down (or vice versa) as the external 
field is changed. This jump can trigger instabilities in 
other pillars, forming an avalanche. We believe that in 
order for multicycles to be seen, it is essential for the 
interaction between pillars to be sufficiently strong to 
cause avalanches; in the extreme case, when the pillars 
are independent, it is clear that a one-cycle hysteresis 
loop would be seen. However, avalanches are not suffi- 
cient to produce disorder: for disordered nearest neigh- 
bor Ising ferromagnets, the phenomenon of Return Point 
Memory (RPM) can be proved, precluding mul- 

ticycles. This emphasizes the need for frustration. 

This is the first paper, as far as we are aware, that stud- 
ies adiabatic hysteresis loops in magnetic systems using 
the more fundamental LLG equations rather than sim- 



plified relaxational dynamics. Experience from critical 
phenomena might lead one to believe that the difference 
between this approach and previous work would be triv- 
ial. However there is an important physical difference 
that we believe has been overlooked. LLG dynamics de- 
stroy symmetry under global spin flip, even though the 
Hamiltonian is symmetric under this operation. This re- 
sult, which will be discussed further in the next section of 
this paper, has significant experimental implications |l3| . 

In the next section of this paper, the dynamical equa- 
tion used in the numerical simulations is introduced, and 
various terms in the model Hamiltonian are calculated. 
Details of the numerics are given in Section III, and the 
results thereof are presented in Section IV. 



II. CLASSICAL SPIN DYNAMICS AND THE 
MODEL HAMILTONIAN 

Microscopically, the evolution of classical spins is de- 
scribed by the Landau-Lifshitz-Gilbert (LLG) equation 
of motion yjj. The LLG equation is the simplest equa- 
tion describing micromagnetic dynamics which contains 
a reactive term and a dissipative term: 

ds ^ 

— = -7is X S - 72S X (s X B), (1) 

where s is a microscopic spin, B is the local effective 
field, 71 is a precession coefficient, and 72 is a damping 
coefficient. The effective field \s B = — ^ + C> where Ti, 

is the Hamiltonian and C, represents the effect of thermal 
noise. Terms in the Hamiltonian will be discussed and 
computed later in this section. 

The reactive term of the LLG equation describes the 
precession of the spin about its local field, with the an- 
gle between the two remaining constant. (The coefficient 
of the reactive term 71 will be set to unity throughout 
this paper unless otherwise noted.) The dissipative term 
aligns the spin with its local effective field. The cross 
products in the dissipative term ensure that only the 
tangential component of the field causes damping, since 
the length of s cannot change. The relaxation time is 
inversely related to the damping coefficient 72. Reason- 
able approximations for 72 are difficult to obtain, but it 
will be shown that the hysteresis multicycle phenomenon 
studied in this paper is present for a large range of 72 . 

With current technology, nanomagnetic pillars that are 
approximately 50nm wide and lOOnm tall can be made of 
ferromagnetic materials such as nickel j2| • For such small 
pillars, it is found that the ferromagnetic coupling be- 
tween the atoms dominates the antiferromagnetic dipo- 
lar interactions. Thus the entire pillar consists of a sin- 
gle magnetic domain. Using the lattice constant of nickel, 
each pillar holds approximately 10^ atoms, allowing us to 
treat the pillar as a continuous magnetic medium. Edge 
effects such as splaying near the boundaries are neglected, 
and the pillar is treated as a saturated nanomagnet with 
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uniform magnetization. Each single-domain nanomagnet 
can be viewed as a single degree of freedom: a magnetic 
moment of fixed magnitude, whose orientation represents 
the direction of the magnetization The time evolu- 
tion of this magnetic moment has the same structure 
as the micromagnetic LLG equation: a reactive part, 
— 7is X B, and a dissipative part, — 72s x (sx B), although 
72 is different from its microscopic value. (Henceforth s 
will denote a unit vector in the direction of the magnetic 
moment of a pillar, rather than an individual spin.) As 
before, the field B is given by the large number of 

spins evolving in unison in each pillar allows the thermal 
noise ^ to be neglected. 

To complete the specification of the dynamics of the 
magnetic moments through Eq. ^ , the Hamiltonian has 
to be calculated in terms of the magnetic moment of each 
pillar. The various terms in the Hamiltonian are dis- 
cussed in the following paragraphs. 

First, the geometry of the pillars introduces a shape 
anisotropy term in the Hamiltonian: 

HsA^-d^Y^sl,, (2) 

i 

where is a constant to be calculated in the next para- 
graph and Sz^i is the z-component of the magnetic mo- 
ment of the i^^ pillar. Shape anisotropy energy is present 
because of the dipolar interactions between the individual 
spins within a pillar. Qualitatively, if the magnetization 
of a tall skinny pillar is vertical, the spins are predomi- 
nantly lined up "head to toe" . This configuration has a 
lower energy than when the magnetization is horizontal, 
in which case the spins are predominantly side by side. 
For a short wide disk, the effect is clearly reversed. 

For the case of tall pillars, the shape anisotropy reduces 
the magnetic moment to an almost Ising like variable, 
that can (approximately) only point up or down. The 
dynamics of anisotropic and isotropic spins are qualita- 
tively different, with avalanche phenomena more likely to 
occur in the former than the latter. As mentioned earlier, 
we believe that avalanches are necessary for hysteretic 
multicycles. Note that even when the shape anisotropy 
is large, we evolve each magnetic moment according to 
Eq.(^) rather than as an Ising variable, i.e. with an 
orientation that evolves continuously with time, though 
the shape anisotropy causes rapid transitions from up to 
down states. 

Deriving the form of Ti-sA and the value of re- 
quires solving a magnetostatic problem. The energy of 
the field due to the microscopic spins in a single pillar is 
W = 2^ / d^x\B{x)\^ , where B[x) is the magnetic field 
at X due to the spins. Through Ampere's law and vec- 
tor calculus manipulations, we can rewrite this in terms 
of the magnetization, M{x). For uniform magnetization, 
the result can (up to an additive constant) be converted 
to a surface integral, similar to electrostatics, with the 
self-energy of the magnetic surface "charge" to be calcu- 
lated. For cylindrical pillars, which we consider in the 



TABLE I: Calculated shape anisotropy coefficients d^, for pil- 
lars with radius _R=30 nm and different aspect ratios h/R. 
The ratios between crystal anisotropy coefficients and are 
also given. 



h/R 


d, (xlO-i*J) 


\Kx/d.\ 


\K2/dA 


0.5 


-2.866 


0.0740 


0.0296 


1 


-2.689 


0.1577 


0.0631 


2 


0.769 


1.102 


0.441 


3 


5.754 


0.2211 


0.0885 


4 


11.30 


0.1501 


0.0600 


5 


17.11 


0.1240 


0.0496 


10 


47.41 


0.0894 


0.0358 



rest of this paper, if the magnetization has magnitude 
Mo and makes an angle 9 to the vertical, and R and h 
are the radius and height of the pillars, the final result is 

W fioM^R^F{^) cos^ e (3) 
R 

up to an additive constant independent of 6. Compar- 
ing with Eq.Q, we see that d^ = -noM^R^F{h/R), 
where F{h/R) is a function of the aspect ratio that can 
be evaluated numerically. With Mq equal to the satu- 
ration magnetization for nickel, 4.84 x W^A/m Tg'I, the 
values of dz for different sized pillars are given in Tabled 
If the pillars are ellipsoidal instead of cylindrical, d, can 
be obtained analytically instead of numerically [isL Il6| . 

A second form of anisotropy energy is caused by the 
crystal structure of Ni. As mentioned earlier, the crystal 
axes give a preferential direction to the magnetization, 
independent of the shape of the pillar. Nickel has a face- 
centered cubic structure which tends to align spins in 
the [111] direction. If cky,ii and Uza are the direc- 
tion cosines of the magnetization of the i^^ pillar to its 
(x, y, z) crystal axes, the crystal anisotropy energy can be 
expanded in powers of a. The first two terms are |l5lll6| 

i 

where an additive constant has been dropped. The ma- 
terial parameters Ki and K2 can be obtained by multi- 
plying the experimentally obtained energy densities with 
the volume of the pillars. Crystal anisotropy energy 
densities have approximate values of —5 x 10^ J/m^ and 
—2 X 10^ J/m^ for Ki/vol and K2/V0I respectively [Tfij . 
Though the ratio of K2/K1 is ~ 0.4, the first term dom- 
inates, since it has 2 less powers of a. 

In addition to the shape and crystal anisotropies that 
affect each pillar by itself, there is dipolar coupling be- 
tween pillars. Microscopically, this is similar to the shape 
anisotropy energy, except that it arises from interactions 
between spins on different pillars. The resultant interac- 
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tion energy is of the form 

T-Cdip ^ ^ Si ■ A{r,j) ■ Sj . (5) 

Aifij) is a second rank tensor that depends only on the 
separation of the piUars. The elements of A{fij) are de- 
termined by numerically solving integrals similar to the 
integrals for the shape anisotropy energy. 

The last term in the Hamiltonian is due to the ex- 
ternal magnetic field, which we take to be in the z- 



direction. The form of this term is the conventional one, 
T^ext = —Be^'^^Sz.i- Hysteresis occurs as is varied 
adiabatically, with the system evolving according to the 
LLG equation. 

In summary, the full Hamiltonian has four terms: 
shape anisotropy, crystal anisotropy, external field, and 
dipole-dipole interaction. The first two are properties of 
the pillars individually, the external field term is the term 
that is adiabatically changed to observe hysteresis, and 
the dipole term is an interaction between pillars: 



Kl, 4 

— K, 



2 2 
i^y.i z,' 



MA 



(6) 



In the numerics, the actual values of the coefficients in 
the Hamiltonian are inconsequential, and only the ratios 
of terms are relevant. Table HI shows the results of the 
calculations described above for c?^, Ki, and if 2. Ev- 
idently, for the dimensions of the nanomagnetic pillars 
of interest, the shape anisotropy term is larger than the 
crystal anisotropy. The dipolar coupling is also small 
compared to dz for the lattice spacings of interest. In the 
simulations, all energies are normalized to dz. Although 
dz is dominant for pillars with aspect ratios of interest, 
the other terms must be included in the Hamiltonian be- 
cause they affect the dynamics qualitatively. Without 
the dipole term, each pillar would be isolated. The crys- 
tal anisotropy term introduces quenched randomness in 
the system, and determines the order in which the s^'s 
flip when the magnetic field is changed; its importance 
has been discussed towards the end of Section I. Thus the 
Hamiltonian of Eq.® has all the important terms that 
have to be kept. 

As mentioned in Section I, although the Hamiltonian 
of Eq.® is invariant if all the siS are flipped (along with 
the external magnetic field) , the dynamics of Eq. are 
not. In Eq.lQJ, the left hand side and the dissipative 
term on the right hand side change sign, but the reactive 
term does not. Therefore the spin inversion symmetry, 
although relevant to equilibrium static properties, does 
not apply to the non equilibrium dynamics appropriate 
for hysteresis. In particular, the two sides of the major 
hysteresis loop are not complementary to each other. 



III. NUMERICS 

As mentioned earlier, the pillars are modeled as sin- 
gle degrees of freedom which follow the LLG equation of 
motion. The effective field for each magnetic moment in 
the LLG equation is the 'spin' derivative of the Hamil- 
tonian of the previous section. Pillars are placed on a 



2-d triangular lattice, to maximize the frustration of the 
dipolar bonds. All systems studied are 4x4 lattices with 
open boundary conditions. The positions of the pillars 
are ix + j{x /2 + y^/?,/2) with i,j = 0, ...,3. The orien- 
tation of the crystallographic axes is separately and ran- 
domly chosen for each pillar. Depending on the choice of 
these random orientations, a sample can have multicy- 
cles of various orders m, or a simple hysteresis loop (i.e. 
m = 1). Square lattices are also considered. 

The dimensions of the cylinders and the separation be- 
tween them, the external field range, and the damping 
coefficient 72 are input parameters. These are used to 
calculate dz, Ki, K2, and the elements of A{rij). All the 
input parameters can be adjusted to maximize the oc- 
currence of multicycles, except 72. Since 72 is a property 
of the material, but is unknown, we make sure that the 
results reported here are valid over a wide range of 72 . 

Numerical modeling of the adiabatic field variation is 
straightforward. The external field is lowered or raised by 
a small field step, SB^- To optimize speed, the field step 
SBf, is adjusted adaptively, since a small step is required 
during avalanches. The effective field is then calculated 
and the system evolves by a small time step, St, with this 
field. This time evolution is repeated, without changing 
Be- Numerical integration of the LLG equation is imple- 
mented using the fourth order Runge-Kutta algorithm. 
Once the system "settles" to a stationary state, the ex- 
ternal field is changed again. Waiting for the system to 
reach a stationary state is equivalent to varying Be slower 
than all the dynamics of the system, i. e. adiabatically. 

The requirement for settling is that the configuration 
after the evolution by a time step St is essentially the 
same as the configuration before. In practice, some nu- 
merical tolerance is allowed, and the initial and final con- 
figurations must differ by less than this tolerance. The 
sums i -I- Ssy i + Ssl i for each i and J2i Ssza must all 
be less than I0~^^ in one time step for the system to be 
considered stationary. The results reported here are in- 
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FIG. 1: Major hysteresis loop (solid curve) for a 4 x 4 tri- 
angular lattice of pillars. The steps demonstrate avalanching 
dynamics. In order to see that the two halves of the loop are 
not complementary, the dashed curve shows the same hystere- 
sis loop, with M —f —M and B — > —B; the solid and dashed 
curves clearly do not coincide. 



FIG. 2: A hysteresis 2-cycle, starting at Be = -BT"'^ . The 
solid curve is a hysteresis loop after one cycle of the external 
field. The dashed curve is the hysteresis loop after the sec- 
ond cycle. Another sweep of the external field would retrace 
the solid curve, indicating that this particular realization of 
randomness undergoes a 2-cycle. B™"^ is less than the sat- 
uration field Bl"*. (If B™"^ were increased beyond SI"*, no 
multicycles would be found.) 



sensitive to a reduction of SB^, St, or the tolerance, and 
therefore represent adiabatic field variation with contin- 
uous time dynamics. 

Starting from a large positive (so that aU the pillars 
are magnetized upwards) the external field is lowered and 
cycled adiabatically over a range [— i?™"^]. The 
configuration {si} is compared at iJ™"^ after each cy- 
cle. If {si} is the same after every m occurrences of 
Be — i?™"^, the system is in an m-cycle. Similar to the 
condition for settling, the configurations match up to a 
tolerance; we have verified in numerous cases that the 
tolerance does not introduce spurious multicycles. A tol- 
erance of 10""* for each component of the magnetization 
was found to be sufficient. Initially, the system undergoes 
a transient period of a few cycles of B^ before reaching a 
limit cycle. 

Figure ^ shows the major hysteresis loop for a sample 
realization of randomness. Since all the pillars are mag- 
netized vertically each time Bg = ±5™"^ = ±oo, m is 
trivially equal to 1. However, one can see the avalanch- 
ing dynamics characteristic of this system, and the fact 
— discussed earlier — that the two halves of the major 
loop are not complementary. Figure |21 shows a hysteresis 
2-cycle. 

IV. RESULTS 

Using an algorithm that performs the operations of 
the previous section, we search through a large number 
of realizations of randomness to find regions in parameter 
space where the probability of finding multicycles is high. 
The parameters in the model are radius R and height h of 
the pillars, _B™°^, and the damping coefficient 72. If the 
lattice spacing, R and h are all scaled by a factor A, all 
the terms in the Hamiltonian are scaled by , which does 



not affect the ratios of the terms. Accordingly, the lat- 
tice spacing can be set to 100 nm without loss of general- 
ity. Given a set of parameters, the algorithm determines 
the periodicity m for a given realization of randomness. 
By classifying the periodicity for a large number of re- 
alizations, we obtained the approximate probability for 
finding an m-cycle as a function of to. 

The easiest parameter to vary experimentally is iJ™"^ . 
As mentioned in the previous section, when is 
too large or too small, multicycles will not be present. 
We find that multicycles can be roughly optimally found 
when 5™°^ approaches the saturation field but not 
greater. is different for every realization of random- 
ness, therefore the optimal field can only be determined 
by scanning over various values of B™"^. The range of 
j^max -^yjjgj-g occurrence of multicycles is apprecia- 
ble depends on the pillar dimensions. For i?=30 nm and 
h—l^Q nm, the probability of finding multicycles when 
j^max optimal is ~ 2 — 3 times the probability when 
j^max percent from the optimal field. In general, 

for systems in which the probability of finding multicy- 
cles is small, the range of iJ™<*^ is narrow. The narrow 
range in B™""^ is not an obstacle to finding multicycles, 
due to the ease of tuning this parameter experimentally. 

The damping coefficient of the LLG equation, 72, can- 
not be easily calculated. In fact, different experimental 
environments could allow for a large range of 72 . Because 
of our inability to obtain a reasonable and realistic ap- 
proximation for 72 , we run searches for a wide range of 72 
(relative to 71). The results show that multicycles exist 
for very small 72 to essentially infinite damping. (The 
large damping limit is implemented by setting 71 =0 and 
keeping 72 finite.) Small values of 72 tend to give more 
multicycles, as one might expect: the probability of find- 
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TABLE II: Approximate probabilities of finding an m-cycle TABLE III: Approximate probabilities of finding an m-cycle 
for systems with pillars of difi'erent radii. The lattice spacing for systems with pillars of difi'erent aspect ratios h/ R. The 
is 100 nm and the aspect ratio, h/ R, is 5. lattice spacing is 100 nm and the radius is fixed at 30 nm. 



R (nm) 


P™.=2 


Pm=3 


-Pm— 4 


Pm>A 


Pm>l 


h/R 


-Pm=2 




-Pm— 4 






10 

















.5 

















20 


.08 


.02 








.1 


1 

















25 


.14 


.12 


.02 


.04 


.32 


3 


.16 


.04 


.02 


.12 


.34 


30 


.2 


.06 


.04 


.16 


.46 


4 


.12 


.08 


.02 


.16 


.38 


35 


.22 


.22 


.04 


.12 


.6 


5 


.2 


.06 


.04 


.16 


.46 


40 


.22 


.12 


.1 


.18 


.62 


6 


.16 


.1 


.06 


.16 


.48 


45 


.28 


.12 





.10 


.50 


10 


.24 


.2 


.06 


.1 


.6 



ing multicycles increases by a factor of 1.5 when 72 is 
reduced from 71 to ~ O.I71. Unexpectedly, the multicy- 
cle probability also seems to increase slightly when 72 is 
larger than 71. When the dynamical equation is strictly 
dissipative (71=0), the multicycle probability is compa- 
rable to the probability when 72 ^ 71 • There are no clear 
trends in the distribution of m when 72 is changed. 

Because of the difficulty in calculating 72 , we conser- 
vatively set 72 to a value where the probability of find- 
ing multicycles is approximately minimal. As mentioned 
in the preceding paragraph, this occurs when 72 ~ 7i ■ 
With at its optimal value, and 72 = 71 = 1, prob- 

abilities of finding multicycles for systems with pillars of 
different R and h are found. The different terms in the 
Hamiltonian scale differently as R and h are varied. Ta- 
ble ^1 shows results for systems of the same aspect ratio 
{h/R = 5) and different radii. The maximum radius for 
a pillar is 45 nm, with the lattice spacing 100 nm. From 
Table ^ one can see that systems of pillars with larger 
radii generally display multicycles more often. 

The occurrence of multicycles depends largely on the 
aspect ratio of the pillars. Avalanches tend to occur only 
for pillars where h/R is large; accordingly, systems with 
disk-shaped pillars, i.e. with negative dz, do not display 
multicycles. In fact, even when dz is positive, multicycles 
are only found for a sufficiently large shape anisotropy 
energy. Table ITTll shows numerical results for systems of 
pillars with a 30 nm radius and various aspect ratios. 
Multicycles are only likely to be found for systems with 
long pillars. No multicycles are found for disk-shaped 
pillars as expected. 

Systems with pillar vacancies are also studied. Vacan- 
cies are implemented by introducing a finite probability 
for a pillar to be missing at every lattice site. One might 
expect that these vacancies would introduce more ran- 
domness in the system, thereby increasing the number of 
multicycles. For pillars with R=30nm and h=150nm, the 
probability of finding a multicycle is ~ 0.46 without any 
vacancies. When the probability of having a vacancy at 
a site is small 0.2), the number of multicycles drops 
by about 40 percent. When the vacancy probability in- 
creases to > 0.5, the probability of finding multicycles 
decreases to less than 0.1. This decrease in probability 



could be due to the decrease in number of pillars. We con- 
clude that, contrary to what one might expect, random 
vacancies do not increase the probability of multicycles . 

If a square lattice is used instead of a triangular one, 
the dipolar couplings between nearest neighbor pillars 
are not frustrated. In a checkerboard pattern, all nearest 
neighbor bonds would be satisfied, but the next near- 
est neighbor bonds (and certain neighbors further apart) 
would not be satisfied. We find that the amount of frus- 
tration in the square lattice is sufhcient for multicycles to 
be found, though the probability of finding a multicycle 
is significantly less on a square lattice than on a triangu- 
lar lattice. The number of multicycles found on a square 
lattice is approximately half of the number found on a 
triangular lattice. 

A possible mechanism for increasing frustration and 
disorder is to introduce random ferromagnetic couplings 
between pillars. These bonds could be manufactured by 
building one large pillar instead of two small ones. The 
large pillar would have two domains that are coupled fer- 
romagnetically. By randomly assigning the bonds with 
some probability, more disorder can be introduced in ad- 
dition to that from the crystal anisotropy. We did not 
study how these ferromagnetic bonds affect the number 
of multicycles. 

One interesting question is whether exact Return Point 
Memory survives when it is extended beyond the 
random field Ising model with purely ferromagnetic in- 
teractions to continuous time vector models. To answer 
that, we use the LLG model with nearest neighbor ran- 
dom ferromagnetic interactions and the same crystalline 
and shape anisotropy model that is used above (without 
dipole coupling). We apply a random field and search 
for violations of RPM with different random seeds. We 
start at high fields and go to a minimum field of —1.7 
and record the spin configuration. Then we go up to a 
field 1.7 and back down to —1.7. We find that for a 4 x 4 
square lattice of spins, that of order 1 percent of sys- 
tems violate RPM because the initial and final minimum 
configurations and total magnetization are substantially 
different. Violation of RPM is seen both with and with- 
out the precessional term in Eq.(Q. This shows that it 
is not possible to extend the proof of RPM to continuous 
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vector models. 



V. CONCLUSIONS 

In this paper, we have investigated the feasibihty of 
observing multicycles and non-complementary hystere- 
sis loops in a candidate experimental system: that of 
cylindrical magnetic nano-pillars arranged on a lattice. 
Wc have performed realistic numerical simulations of this 
system, by calculating the magnetic interactions between 
the pillars and then employing continuous spin dynamics 
and the Landau Lifshitz Gilbert equation to obtain their 
time evolution. Using physically appropriate parameters, 
we have shown that there is often multicycle hysteretic 
behavior, i.e. a periodic adiabatic external magnetic field 
causes a subharmonic steady state response in the mag- 
netization. We have also shown that, even though the 



Hamiltonian is invariant under spin and field reversal, 
the dynamics are not, so that the two branches of the 
major hysteresis loop are not complementary. This re- 
sult cannot be obtained with an Ising approximation to 
the system and its dynamics. 

Further implications of the non-complementary nature 
of LLG dynamics will be studied in future research. 

Because systems of this kind are currently the sub- 
ject of much experimental investigation, we believe that 
it would be fruitful to attempt to observe the unusual 
behavior predicted here. 
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